home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / vfft / vfftpass.c < prev   
Encoding:
C/C++ Source or Header  |  1992-06-01  |  10.5 KB  |  411 lines
  1. /*
  2. %    VFFTPASS . C
  3. %
  4. %    Symmetric Filter for VFFT image transform.
  5. %
  6. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  7.  
  8. This software is copyright (C) by the Lawrence Berkeley Laboratory.
  9. Permission is granted to reproduce this software for non-commercial
  10. purposes provided that this notice is left intact.
  11.  
  12. It is acknowledged that the U.S. Government has rights to this software
  13. under Contract DE-AC03-765F00098 between the U.S.  Department of Energy
  14. and the University of California.
  15.  
  16. This software is provided as a professional and academic contribution
  17. for joint exchange. Thus, it is experimental, and is provided ``as is'',
  18. with no warranties of any kind whatsoever, no support, no promise of
  19. updates, or printed documentation. By using this software, you
  20. acknowledge that the Lawrence Berkeley Laboratory and Regents of the
  21. University of California shall have no liability with respect to the
  22. infringement of other copyrights by any part of this software.
  23.  
  24. For further information about this notice, contact William Johnston,
  25. Bld. 50B, Rm. 2239, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  26. (wejohnston@lbl.gov)
  27.  
  28. For further information about this software, contact:
  29.     Jin Guojun
  30.     Bld. 50B, Rm. 2275, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  31.     g_jin@lbl.gov
  32.  
  33. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  34. %
  35. % compile:    cc -O -o DEST/vfftpass vfftpass.c -lscs1 -lccs -lhips -lm
  36. */
  37. char    usage[]="options\n\
  38. [-base #]    base value. [Default=0.0], good max=>2.0\n\
  39. [-c]        convex curve. Default is concave curve\n\
  40. [-hi]        high-pass filter. Default is low-pass\n\
  41. [-f #]        elastic scale factor\n\
  42. [-F #]        number of Frames in sample [16]\n\
  43. [-G]        generate plot data of curves only\n\
  44. [-H #]        high pass adjust\n\
  45. [-L]        Linear filter\n\
  46. [-min #]    minimum value at bottom [0.01]\n\
  47. [-neg]        negate filter (1D -hi)\n\
  48. [-r #]        radius factor. default is 1 (100%). Range 0.1 - 1.0\n\
  49. [-v]        verbose\n\
  50. [<] image [> [-o] VFFT_file]\n";
  51. /*
  52. %
  53. % AUTHOR:    Guojun Jin - Lawrence Berkeley Laboratory    5/1/91
  54. */
  55.  
  56. #include "vfft_fil.h"
  57.  
  58. U_IMAGE    uimg;
  59.  
  60. bool    dimens, curve_dir, hi_pass, Table, Msg;
  61. VType    *itemp;
  62. Filter    *lkt[3], flr=0.01, sc;
  63. int    dimen1len, h_rows, h_frm, sframes=16,
  64.     r_radius, c_radius, f_radius;
  65. MType    fsize, vsize;
  66.  
  67.  
  68. main (argc, argv)
  69. int    argc;
  70. char**    argv;
  71. {
  72. bool    dflag=0, neg=0, vflag, sample=0;
  73. int    f, j;
  74. Filter    base=0, Hbase=.005;
  75. Filter    radius_rate=1., *filter3d, *filter_work;
  76. register int    i;
  77.  
  78. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, -1, *argv, "S20-1");
  79.  
  80. for (j=1; j<argc; j++)
  81.     if (*argv[j] == '-'){
  82.     f=1;
  83.     switch (argv[j][f++]) {
  84.     case 'h':
  85.         hi_pass++;    break;
  86.     case 'L':
  87.         curve_dir=CURVE_LINEAR;    break;
  88.     case 'c':
  89.         curve_dir=CURVE_CONVEX;    break;
  90.     case 'b':
  91.         base = GValue(0);    break;
  92.     case 'f':
  93.         sc = GValue(0);        break;
  94.     case 'F':
  95.         if (GValue(0))
  96.             sframes = atoi(argv[j]+f);
  97.         break;
  98.     case 'G':
  99.         Table++;    break;
  100.     case 'H':
  101.         if (GValue(0))
  102.             Hbase = atof(argv[j]+f);
  103.         break;
  104.     case 'm':
  105.         if (GValue(0))
  106.             flr = atof(argv[j]+f);
  107.         break;
  108.     case 'n':
  109.         neg++;    break;
  110.     case 'r':
  111.         if (GValue(0))
  112.             radius_rate = atof(argv[j]+f);
  113.         if (radius_rate < 0 || radius_rate > 1.0)
  114.             radius_rate = 1.0;
  115.         break;
  116.     case 'v':    Msg++;
  117.         break;
  118. #ifdef    IBMPC
  119.     case 'i':if (GValue(1) && !(in_fp=freopen(argv[j]+f, "rb", stdin)))
  120.             goto    ferr;
  121.         break;
  122. #endif
  123.     case 'o':if (GValue(1) && freopen(argv[j]+f, "wb", stdout))    break;
  124. ferr:        message("File %s can't be opened", argv[j]);
  125.     default:
  126.         usage_n_options(usage, j, argv[j]);
  127.     }
  128.     }
  129.     else if ((in_fp=freopen(argv[j], "r", stdin)) == NULL)
  130.     syserr("can't open frame file %s", argv[j]);
  131.  
  132. io_test(fileno(in_fp), sample=!iset);
  133.  
  134. if (sample) {    /* no input file */
  135.     uimg.width = uimg.height = 128;
  136.     uimg.frames = sframes;
  137.     uimg.pxl_in = 8;
  138.     uimg.in_form = IFMT_VFFT3D;
  139.     (*uimg.std_swif)(FI_INIT_NAME, &uimg, "VPASS", 0);
  140. }
  141. else    {
  142.     (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  143.     if (uimg.in_form < IFMT_VFFT3D || uimg.in_form  > IFMT_VVFFT3D)
  144.         syserr("image is not a VFFT format");
  145.     vflag = uimg.in_form==IFMT_VVFFT3D || uimg.in_form==IFMT_DVVFFT3D;
  146. }
  147. uimg.o_form = uimg.in_form;
  148. uimg.pxl_out = uimg.pxl_in;
  149.  
  150. if (!vflag)    /* the VVFFT must be in 3D */
  151.     dimens = !(uimg.in_form & 1); /* if it's 2D */
  152.  
  153. dimen1len = (cln>>1) + 1;
  154. h_rows = row >> 1;
  155. h_frm = frm >> 1;
  156. fsize = row*cln;
  157. vsize = row*dimen1len;
  158. f_radius = radius_rate * h_frm;
  159. r_radius = radius_rate * h_rows;
  160. c_radius = radius_rate * dimen1len;
  161. if (f_radius < 1)    f_radius = 1;
  162. if (r_radius < 1)    r_radius = 1;
  163. if (c_radius < 1)    c_radius = 1;
  164.  
  165. lkt[lktCol] = zalloc(cln, sizeof(Filter), "lktcol");
  166. lkt[lktRow] = zalloc(row, sizeof(Filter), "lktrow");
  167. lkt[lktFrm] = zalloc(frm, sizeof(Filter), "lktfrm");
  168.  
  169. if (hi_pass)    sc = -sc;
  170.  
  171. new_curve(lkt[lktCol], c_radius, neg, curve_dir, "samples");
  172. new_curve(lkt[lktRow], r_radius, neg, curve_dir, "lines");
  173. new_curve(lkt[lktFrm], f_radius, neg, curve_dir, "layers");
  174. if (hi_pass){
  175. register Filter    *lktp=lkt[lktFrm];
  176.     for (i=f_radius; i--;)
  177.         lktp[i] = 1. - lktp[i] + .01;
  178.     lktp[0] = Hbase;
  179. /*    for (i=r_radius, lktp=lkt[lktRow]; i--;)
  180.         lktp[i] = 1. - lktp[i];
  181.     for (i=c_radius, lktp=lkt[lktCol]; i--;)
  182.         lktp[i] = 1. - lktp[i];
  183. */
  184.     if (Msg)
  185.         dump_lkt(lktp, f_radius),
  186.         msg("\n%d hi_pass => floor=%f\n", f_radius,flr);
  187. }
  188. if (Table)    exit(0);    /* only send plot data to stdout */
  189.  
  190. filter3d = nzalloc(sizeof(*filter3d)<<1, r_radius*c_radius, "filter3d");
  191. filter_work = filter3d + r_radius * c_radius;
  192. ibuf = nzalloc(uimg.pxl_in, vsize, "ibuf");
  193. obuf = zalloc(uimg.pxl_in, vsize, "obuf");    /* be zeroed */
  194.  
  195. {
  196. register Filter    *filterp=filter3d, *llktp=lkt[lktCol];
  197.     memcpy(filterp, llktp, c_radius*sizeof(*filterp));
  198.     filterp += c_radius;
  199.     for (i=1; i<r_radius; i++){
  200.     register Filter    scale = lkt[lktRow][i];
  201.         for (j=0; j<c_radius; j++)
  202.         *filterp++ = scale * llktp[j];
  203.     }
  204.     if (hi_pass){    /* reverse low to high in 2D instead of in 1D    */
  205.         filterp = filter3d;
  206.         for (i=r_radius*c_radius; i--;)
  207.             filterp[i] = 1. - filterp[i];
  208.         filterp[0] += Hbase;
  209.     }
  210. }
  211. (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  212.  
  213. if (sample){
  214. register Filter    *filterp=filter_work;
  215.  
  216.     {
  217.     register float    *ibp=ibuf;
  218.         for (i=vsize << 1; i--;)
  219.         ibp[i] = 256.;
  220.     }
  221.     for (f=h_frm; f--;){
  222.     register Filter    scale = lkt[lktFrm][f];
  223.         for (i=r_radius*c_radius; i--;)
  224.             filterp[i] = filter3d[i] * scale;
  225.         filtering(ibuf, obuf, base, filterp, f+1);
  226.     }
  227.     if (frm & 1)
  228.         filtering(ibuf, obuf, base, filter3d, f+1);
  229.     for (f=0; f<h_frm; f++){
  230.     register Filter    scale = lkt[lktFrm][f];
  231.         for (i=r_radius*c_radius; i--;)
  232.             filterp[i] = filter3d[i] * scale;
  233.         filtering(ibuf, obuf, base, filterp, f+h_frm);
  234.     }
  235.     exit(0);
  236. }
  237.  
  238. if (dimens){
  239.     message("%s: 2 dimension filter - %d frames\n", Progname, frm);
  240.     for (f=0; f<frm; f++){
  241.         i = upread(ibuf, uimg.pxl_in, vsize, in_fp);
  242.         if (i != vsize)
  243.             syserr("f%d [%d] read %d", f+1, vsize, i);
  244.         filtering(ibuf, obuf, base, filter3d, f+1);
  245.     }
  246. }
  247. else    {
  248. register Filter    *filterp=filter_work;
  249.  
  250.     if (hi_pass){
  251.     register Filter    scale = lkt[lktFrm][0];
  252.         for (i=r_radius*c_radius; i--;)
  253.             filterp[i] = filter3d[i] * scale;
  254.     }
  255.     else    filterp = filter3d;
  256.     i = upread(ibuf, uimg.pxl_in, vsize, in_fp);
  257.     if (i != vsize)
  258.         syserr("f%d [%d] read %d", 1, vsize, i);
  259.     filtering(ibuf, obuf, base, filterp, 1);
  260.  
  261.     for (f=1; f<f_radius; f++){
  262.     register Filter    scale = lkt[lktFrm][f];
  263.         filterp = filter_work;
  264.         for (i=r_radius*c_radius; i--;)
  265.             filterp[i] = filter3d[i] * scale;
  266.         i = upread(ibuf, uimg.pxl_in, vsize, in_fp);
  267.         if (i != vsize)
  268.             syserr("f%d [%d] read %d", f+1, vsize, i);
  269.         filtering(ibuf, obuf, base, filterp, f+1);
  270.     }
  271.  
  272.     for (; frm>1 && f<frm-(f_radius<<1); f++){    /* maybe send all 0 frames */
  273.         i = upread(ibuf, uimg.pxl_in, vsize, in_fp);
  274.         if (i != vsize)
  275.             syserr("f%d [%d] read %d", f+1, vsize, i);
  276.         filtering(ibuf, obuf, base, filterp, f+1);
  277.     }
  278.  
  279.     for (; f<frm-1; f++){
  280.     register Filter    scale = lkt[lktFrm][frm-1-f];
  281.         filterp = filter_work;
  282.         for (i=r_radius*c_radius; i--;)
  283.             filterp[i] = filter3d[i] * scale;
  284.         i = upread(ibuf, uimg.pxl_in, vsize, in_fp);
  285.         if (i != vsize)
  286.             syserr("f%d [%d] read %d", f+1, vsize, i);
  287.         filtering(ibuf, obuf, base, filterp, f+1);
  288.     }
  289.  
  290.     i = upread(ibuf, uimg.pxl_in, vsize, in_fp);
  291.     if (i != vsize)
  292.         syserr("f%d [%d] read %d", f+1, vsize, i);
  293.     if (hi_pass){
  294.     register Filter    scale = lkt[lktFrm][0];
  295.         for (i=r_radius*c_radius; i--;)
  296.             filterp[i] = filter3d[i] * scale;
  297.     }
  298.     else    filterp = filter3d;
  299.     filtering(ibuf, obuf, base, filterp, f+1);
  300. }
  301. }
  302.  
  303. new_curve(lkt, maxdiff, neg, curve, name)
  304. float*    lkt;
  305. int    maxdiff;
  306. char    *name;
  307. {
  308. Filter    rel_val, maxout=1.;
  309. register Filter    scale, vsc=.01*sc/maxdiff + 1;
  310. register int    i;
  311.  
  312. if (maxdiff==1){
  313.     lkt[0] = 1.;
  314.     return;
  315. }
  316. if (curve==CURVE_LINEAR){
  317. register int    rev;
  318.     scale = (maxout-flr)/maxdiff;
  319.     if (neg) for (i=maxdiff; i--;)
  320.         lkt[i] = i * scale + flr;
  321.     else for (i=maxdiff, rev=i-1; i--;)
  322.         lkt[rev-i] = i * scale + flr;
  323. }
  324. else{
  325.     {
  326.     register double    tmp;
  327.         if (vsc != 1.){
  328.         for (i=tmp=maxdiff; i--;)    tmp = tmp*vsc + i;
  329.         scale = (maxout-flr) / tmp;
  330.         }
  331.     else    scale = (maxout-flr) / ((maxdiff-1)*maxdiff >> 1);
  332.     }
  333.     if (!neg)
  334.        if (!curve)
  335.         for(i=maxdiff, rel_val=flr; i--;)
  336.         {
  337.         rel_val += (maxdiff-i) * (scale*=vsc);
  338.         if (rel_val > maxout)    lkt[i]=maxout;
  339.         else    lkt[i] = rel_val;
  340.         }
  341.        else for (i=0, rel_val=maxout; i<maxdiff; i++)
  342.         {
  343.         lkt[i] = rel_val;
  344.         rel_val -= i * (scale*=vsc);
  345.         if (rel_val<flr && !uimg.in_form)    rel_val=flr;
  346.         }
  347.     else if (!curve)
  348.         for(i=0, rel_val=flr; i<maxdiff; i++)
  349.         {
  350.         rel_val += i * (scale*=vsc);
  351.         lkt[i] = rel_val;
  352.         }
  353.        else for (i=maxdiff, rel_val=maxout; i--;)
  354.         {
  355.         lkt[i] = rel_val;
  356.         rel_val -= (maxdiff-i) * (scale*=vsc);
  357.         }
  358.     }
  359. if (Table)    GTable(lkt, maxdiff);
  360. if (Msg)
  361.     dump_lkt(lkt, maxdiff),
  362.     msg("\n%d %s => sc=%f, scale=%f\n", maxdiff, name, sc, scale);
  363. }
  364.  
  365. dump_lkt(lktp, n)
  366. register float    *lktp;
  367. {
  368. register int    i;
  369. for (i=0; i<n;)
  370.    {    message("%10.3f", lktp[i]);    if (!(++i & 7))    mesg("\n");    }
  371. }
  372.  
  373. GTable(lp, max)
  374. register float *lp;
  375. {
  376. register int i;
  377. for (i=0; i<max; i++)
  378.     printf("%d %f\n", hi_pass ? max-i : i, *lp++);
  379. printf("\n\n");
  380. }
  381.  
  382. filtering(i_buf, o_buf, base_v, filterp, f)
  383. VType    *i_buf, *o_buf;
  384. Filter    base_v;
  385. register Filter    *filterp;
  386. {
  387. register float    *ibp=i_buf, *obp=o_buf;
  388. register int    i, j;
  389.  
  390.     if (base_v)
  391.         for (i=r_radius*c_radius; i--;)
  392.         filterp[i] += base_v;
  393.     ibp += hi_pass * ((h_rows-r_radius+1)*dimen1len - c_radius);
  394.  
  395.     for (i=0; i<r_radius; i++){
  396.     register int    dis = (row - 1 - (i<<1)) * dimen1len << 1;
  397.         for (j=0; j<c_radius; j++){
  398.         obp[dis] = ibp[dis] * *filterp;
  399.         *obp++ = *ibp++ * *filterp;
  400.         obp[dis] = ibp[dis] * *filterp;
  401.         *obp++ = *ibp++ * *filterp++;
  402.         }
  403.         obp += dimen1len - c_radius << 1;
  404.         ibp += dimen1len - c_radius << 1;
  405.     }
  406.  
  407.     i = fwrite(o_buf, uimg.pxl_in, vsize, out_fp);
  408.     if (i != vsize)
  409.         syserr("f%d [%d] write %d", f, vsize, i);
  410. }
  411.